#include <stdio.h>
#include <gmp.h>
#include "borwein.h"

#define IT_N_SIZE 32

/* variaveis */

mpf_t npi, old_pi, ak, yk, a_tmp, y_tmp;
mpf_t alpha, beta, gama, tmp1, tmp2, err;
unsigned long int pow2;

void set_init_vals(){
	pow2 = 1;
	
	//a0 = 6-4sqrt(2)
	//y0 = sqrt(2)-1

	mpf_set_ui(yk, 2);
	mpf_sqrt(yk,yk);
	mpf_mul_ui(ak, yk, 4);

	mpf_sub_ui(yk, yk, 1);
	mpf_ui_sub(ak, 6, ak);
}


void calc_yk(){

	//alpha <= (1 - yk^4)^(1/4)
	mpf_pow_ui(alpha, yk, 4);
	mpf_ui_sub(alpha, 1, alpha);
	mpf_sqrt(alpha,alpha);
	mpf_sqrt(alpha,alpha);
	
	//beta <= 1 - alpha
	mpf_ui_sub(beta, 1, alpha);

	//gama <= 1 + alpha
	mpf_add_ui(gama, alpha, 1);

	//y_tmp <= beta/gama;
	mpf_div(y_tmp, beta, gama);
		
}

void update_yk(){
	mpf_swap(yk, y_tmp);
}

/*
ak_pre_calc
	tmp1 <= 1 + yk
	tmp2 <= yk^2
	tmp2 <= tmp1 + y2
	tmp2 <= tmp2 * yk

(yk consumido)
*/

void ak_pre_calc(){
	mpf_add_ui(tmp1, yk, 1); //tmp1 = yk + 1
	mpf_pow_ui(tmp2, yk, 2);  //tmp2 = yk^2
	mpf_add(tmp2, tmp2, tmp1); //tmp2 = 1 + yk + yk^2
	mpf_mul(tmp2, tmp2, yk); //tmp2 = yk(1 + yk + yk^2)
}

/*
ak_end_calc
	tmp1 <= tmp1^4
	tmp1 <= tmp1*ak
	pow2 <= pow2 + 2
	tmp2 <= tmp2*2^(pow2)
	ak2 <= tmp1 - tmp2;
*/


void ak_end_calc(){

	mpf_pow_ui(tmp1, tmp1, 4); //tm
	mpf_mul(tmp1, tmp1, ak);
	pow2 += 2;
	mpf_mul_2exp(tmp2, tmp2, pow2);
	mpf_sub(a_tmp, tmp1, tmp2);
}

void update_ak(){
	mpf_swap(ak, a_tmp);
}

void calc_pi(){
	mpf_swap(npi, old_pi);
	mpf_ui_div(npi,1,ak);
}

int check_erro(){
	mpf_sub(old_pi, npi, old_pi);
    mpf_abs(old_pi, old_pi);
	return (mpf_cmp(err, old_pi) < 0);
}

void save_it(char* name){
	static int it = 0;
	FILE *fd;
	char it_name[IT_N_SIZE];

	sprintf(it_name, "%s-it%d.txt", name, ++it);
    fd = fopen(it_name, "w");
    gmp_fprintf(fd, "%.10000000Ff\n", npi);
    fclose(fd);
}

void init_vars(){
	mpf_init(npi);
	mpf_init(old_pi);
	mpf_init(ak);
	mpf_init(yk);
	mpf_init(a_tmp);
	mpf_init(y_tmp);
	mpf_init(alpha);
	mpf_init(beta);
	mpf_init(gama);
	mpf_init(tmp1);
	mpf_init(tmp2);
	mpf_init_set_str(err, "1e-10000000", 10);
}

void clear_vars(){
	mpf_clear(npi);
	mpf_clear(old_pi);
	mpf_clear(ak);
	mpf_clear(yk);
	mpf_clear(a_tmp);
	mpf_clear(y_tmp);
	mpf_clear(alpha);
	mpf_clear(beta);
	mpf_clear(gama);
	mpf_clear(tmp1);
	mpf_clear(tmp2);
	mpf_clear(err);
}
